On the forces acting on a small particle in an acoustical field in a viscous fluid 



O 
(N 

> 

o 
O 

i> 

(N 
"a 

i 

q 

O 

^ ; 

Oh 



> 

m 
o 



Mikkel Settnes and Henrik Bruus 
Department of Micro- and Nanotechnology, Technical University of Denmark, 
DTU Nanotech Building 345 B, DK-2800 Kongens Lyngby, Denmark 
(Dated: 27 October 2011) 

We calculate the acoustic radiation force from an ultrasound wave on a compressible, spherical 
particle suspended in a viscous fluid. Using Prandtl-Schlichting boundary-layer theory, we include 
the kinematic viscosity of the solvent and derive an analytical expression for the resulting radiation 
force, which is valid for any particle radius and boundary-layer thickness provided that both of these 
length scales are much smaller than the wavelength of the ultrasound wave (mm in water at MHz 
frequencies). The acoustophoretic response of suspended microparticles is predicted and analyzed 
using parameter values typically employed in microchannel acoustophoresis. 



I. INTRODUCTION 

Particles suspended in acoustic fields are subject to 
time-averaged forces from scattering of the acoustic 
waves. Theoretical studies of these forces, knovifu as 
acoustic radiation forces, date back to King in 1934, who 
considered incompressible particles suspended freely in 
an inviscid fluid [1]. In 1955 Yosioka and Kawasima 
extended the analysis to include the compressibility of 
the suspended particles [2]. These results were summa- 
rized and generalized by a simple and physically intuitive 
method by Gorkov in 1962 [3], but limited to inviscid fiu- 
ids and particles smaller than the acoustic wavelength A. 

With recent developments in microfabrication tech- 
nologies allowing for integration of ultrasound resonators 
in lab-on-a-chip systems, the acoustic radiation force has 
received renewed attention as a label- and contact-free 
way to manipulate particles. Several biotechnological ap- 
plications of particle trapping and separation have been 
reported where ultrasound resonances in microchannels 
were used to create acoustic fields giving rise to acous- 
tic radiation forces on suspended particles. Examples 
are on-chip acoustophoretic cell separation devices [4, 5], 
cell trapping [6-8], plasmapheresis [9], forensic analysis 
[10], food analysis [11], cell sorting using surface acoustic 
waves [12], cell synchronization [13], and cell differentia- 
tion [14]. At the same time, substantial advancements 
in understanding the fundamental physics of biochip 
acoustophoresis have been achieved through full-chip 
imaging of acoustic resonances [15], surface acoustic wave 
generation of standing waves [16], multi-resonance chips 
[17], advanced frequency control [18, 19], on-chip inte- 
gration with magnetic separators [20] , acoustics-assisted 
microgrippers [21], in-situ force calibration [22], and au- 
tomated micro-PIV systems [23]. 

Traditionally, the acoustic radiation force has been 
modeled using the inviscid theory of the acoustic radi- 
ation force. This approach is approximately correct for 
particles significantly larger than the thickness 6 of the 
acoustic boundary layer, in which viscosity do play a 
dominant role. For a fluid with kinematic viscosity (or 
momentum diffusivity) v and with an acoustic field of 
angular frequency w, the boundary layer thickness (or 



viscous penetration depth) is the momentum diffusion 
length given by [24, 25] 



0.6 pm. 



(1) 



where the value is for 1-MHz ultrasound in water at room 
temperature. It is therefore expected that particles or 
cells with a radius larger than 3 jim can be described 
fairly accurately by the inviscid theory. However, as the 
technological development pushes for higher accuracy, 
more refined applications, and the handling of smaller 
particles, it becomes relevant to calculate the effects on 
acoustophoresis from viscosity of the solvent. 

We are not the first to analyze how the acoustic ra- 
diation force depends on 5. However, the earlier works 
by Doinikov [26] and by Danilov and Mironov [27] fo- 
cus mainly on developing general theoretical schemes for 
particles of radius a smaller than the wavelength A and 
only provide analytical expressions in the special limits 
S <^ a <^ X and a <C 5 <C A. However, given the length 
scale above, the range of applicability of the published 
expressions for viscous corrections is a priori severely 
limited. The aim of this paper is to provide an analyti- 
cal expression for the viscous corrections to the acoustic 
radiation force on small suspended particles 6^a <^ A, 
and to analyze its implications for experimentally rele- 
vant parameters central for current studies in the field of 
microchannel acoustophoresis of compressible particles in 
liquids. 

We begin by establishing the governing equations for 
acoustophoresis in the framework of second-order pertur- 
bation theory of the Navier-Stokes equation in the acous- 
tic field. Then, following the inviscid analysis by Gorkov 
[3] , we express the radiation force on a particle in terms 
of the far-field solution of inviscid acoustic wave scatter- 
ing theory, extend this solution to the near-field region 
close to the particle, and match it with the solution to 
the incompressible viscous flow problem in the acoustic 
boundary layer of the particle. From this we obtain the 
analytical expression for the acoustic radiation force in 
the viscous case. Finally, we analyze the predictions of 
the theory for experimentally relevant parameter values. 
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II. PERTURBATION EXPANSION OF THE 
GOVERNING EQUATIONS 

The formulation of the governing equations for acous- 
tics in perturbation theory is well known, and the reader 
is referred to the textbooks by Lighthill [28], Pierce [29], 
and Landau & Lifshitz [25]. Briefly and to establish no- 
tation [30], for a given fluid in the absence of external 
forces and for isothermal conditions, the theory is based 
on a combination of the thermodynamic equation of state 
expressing pressure p in terms of density p, the kinematic 
continuity equation for p, and the dynamic Navier-Stokes 
equation for the velocity field v, 



P = P{p) 

dtP = -V-{pv), (2b) 
p^^v = -Vp- p{v-V)v + ryV^t! + l3r] V{V-v), (2c) 



(2a) 
(2b) 



where r] is the dynamic viscosity of the fluid and /3 the 
viscosity ratio typically of the order of unity. Thermal 
effects arc neglected, because the thermal diffusion length 
in liquids is much smaller than the momentum diffusion 
length (or viscous penetration depth) S, see Ref. 26. 

Wc consider a quiescent liquid, which before the pres- 
ence of any acoustic wave has constant density Pq and 
pressure pg- Let an acoustic wave constitute tiny per- 
turbations to first and second order (subscript 1 and 2, 
respectively) in density p, pressure p, and velocity v, 



P 

P ■ 



■ PQ 

■Pa 



Pi +P2> 

Co Pi 



(3a) 
(3b) 
(3c) 



Here, we have introduced the speed of sound Cg of the 
fluid, the square of which is given by the (isentropic) 
derivative Cq = {dp/dp)s, ensuring the useful identity 



Pi = Co Pi, (4) 
and an explicit expression for the compressibility Kg, 
1 dV 1 dp 1 



V dp Pq dp 



(5) 



The first-order perturbation (or linearization) of the 
continuity and Navier-Stokes equation is. 



9tPi = -Po'^-Vi, 



(6a) 

/3?7 V(V-t;i), (6b) 



The first-order acoustic wave equation for p^ is obtained 
by taking the time derivative ^^ of Eq. (6a) followed by 
insertion of Eq. (6b) in the resulting expression. 



d^i 



1 + 



5, 



(7) 



For acoustics fields in the bulk, several times 5 from rigid 
boundaries, the viscous dissipation is negligible because 



of the minute damping coefficient, rjco/ (p^Cq) ^ 1, where 
cj is a characteristic angular frequency of the system. For 
distances within a few times (5 of a rigid wall, the no- 
slip boundary condition forces the velocity of the fluid 
to equal that of the wall, and large velocity gradients 
may occur in Eq. (6b), such that 'qv-^/d'^ > pg^t^i' 
viscosity cannot be neglected. 

The acoustic radiation force is a time-average effect 
that docs not resolve the oscillatory behavior of the 
acoustic fields, so in this work we do not need the full 
second-order perturbation of the governing equations, 
but only their time-average. We assume that after the 
vanishing of transients, any first-order field f{r,t) has a 
harmonic time dependence, 



f{r,t)=fir) e- 



(8) 



and define the time average (X) over a full oscillation 
period r of a quantity X{t) as 



{X) = - f dtX{t). 
Jo 



(9) 



From this we obtain the time-averaged, second-order per- 
turbation of Eqs. (2b) and (2c) in the form 

Po^-(v2) = -'^-(piVi), (10a) 

- V(P2) + ^^(^^2) + /3r?V(V- (t;2)) 

= (Pi^t^^i) +Po((^l•'^K)• 
(10b) 

We note that the physical, real- valued time average {f g) 
of two harmonically varying fields / and g with the com- 
plex representation Eq. (8), is given by the real-part rule 



(/5> = ^ Re 



f{r)g*{r) 



(11) 



where the asterisk denotes complex conjugation. 

Clearly, the time-averaged, second-order fields will in 
general be non-zero, as the non- vanishing time-averaged 
products of first-order terms act as source terms on the 
right-hand side in the governing equations. 

In the inviscid bulk, the first-order flow v-^ is a poten- 
tial flow, see Eq. (6b) with 77 = 0, which when used in 
Eq. (10b) and combined with Eq. (5) leads to 



{P2) = li^oiPi) - Ipoi-"!)- 



(12) 



III. THE ACOUSTIC RADIATION FORCE 

We are analyzing the acoustic radiation force on a com- 
pressible, spherical, micrometer-sized particle of radius 
a suspended in a viscous fluid in an ultrasound field of 
wavelength A (= 1 mm in water at room temperature 
for ui/2tt = 1.5 MHz), thus a <C A. In terms of acoustic 
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(a) _<^,, Far field (b) 




FIG. 1. (a) Sketch of the far-field region r A of an incoming acoustic wave (;/);„ (blue lines) of wavelength A scattering ofi^ a 
small particle (black dot) with radius a <C A, leading to the outgoing scattered wave (j)^^ (red circles and arrows). The resulting 
first-order wave is tpi = -I- (ji^c- (b) Sketch of a compressible spherical particle (yellow) of radius a, compressibility Hp, and 
density pp, surrounded by the incompressible, viscous acoustic boundary layer of width ~ 55 (dark blue) with density pg and 
kinematic viscosity v. Outside is the compressible inviscid bulk (light blue) of compressibility Kq and density Pq. The bulk 
liquid is divided into the near-field region for r <C A, with the instantaneous scattered field (pscit), and the far-field region with 
time-retarded scattered field (jjsdt — t/cq)- 



waves, the microparticle thus acts as a weak point scat- 
terer, which we will treat by first-order scattering theory. 
In response to an incoming wave, described by the oscil- 
latory velocity field v-^^, an outgoing wave v^^ propagates 
away from the particle. For sufiiciently weak waves, the 
first-order acoustic velocity field Vi is given by the sum 

'^i ==^i„ + ^sc- (13) 

Once the first-order scattered field v^^ have been deter- 
mined for a given incoming first-order field v-^^ , the acous- 
tic radiation force f '^'^ on the particle can be calculated 
as the time-averaged second-order forces acting on a fixed 
surface dfl in the inviscid bulk, encompassing the parti- 
cle [3] . Momentum conservation and zero bulk forces en- 
sures that any fixed surface can be chosen. For inviscid 
fluids, F'""^ is the sum of the time-averaged second-order 
pressure (pj) ^^id momentum flux tensor p^i^ViVi'^, 



Tirad 



da 



{(P2)" + Po(("-^i)«i)} (14) 



da 



an 



Po((n-fi)t'i) 



To ease the determination of and Vj^, we use that 
in the inviscid bulk, they can be expressed in terms of 
a velocity potential as Vi — ^(f>i and pi — ^p^d^cfii. 
For a harmonic time dependence, Eq. (6b) implies 



'/'i^-i^Pi, (15) 
and, as sketched in Fig. 1(a), we henceforth write 



V01 



SC 1 

--Vo 



(16a) 
(16b) 
(16c) 



By virtue of Eqs. (7) and (15), 4>i (as well as (f)^^ and 
0j,j.) obeys the inviscid wave equation d^cj) = Cg V^(/) in 
the bulk. As we can use any surface dfl to calculate 
the radiation force F'^^'^, the simplest choice is a sur- 
face in the far-field region r ^ A, where the spherical 
particle of radius a is placed at the center of the coor- 
dinate system, and where r is a position vector. Ac- 
cording to standard scattering theory, the scattered field 
from a point scatterer can be represented by a time- 
retarded multipole expansion. In the far-field region, the 
monopole component and dipole components dominate, 
'?^sc ~ <^mp + "^dp' ^^"i general, these two components 
have the specific forms (t)^^p{r,t) — b{t — r/cg)/r and 
(p^pir, t) = V • \_B{t — r / CQ)/r\, where 6 is a scalar func- 
tion and B a vector function of the retarded argument 
t — r I Cq. In first-order scattering theory, (^^^ must be 
proportional to the first-order fields determined by 
On physical grounds, the only relevant scalar field is the 
density, h ~ pj^^, or equivalently the pressure pj^^, while 
the only relevant vector field is the velocity, B ^ v-^^^. 
Here both p.^^ and v.^^ are evaluated at the particle po- 
sition with time-retarded arguments, and in the far-field 
region must therefore have the form 



0,,(r,t)=-/i 



3po r 



r>A, (17) 



where the particle radius a, the unperturbed density pg, 
and the time derivative are introduced to ensure the 
correct physical dimension of ^g^., namely m^/s. The 
factors 1/3 and 1/2 are inserted for later convenience. 

Before reaching the main goal of the calculation, 
namely the determination of the dimensionless scatter- 
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ing coefficients and the radiation force ir"'a-d jg ex- 
pressed in terms of the incoming acoustic wave (l)-^^ at the 
particle position and the coefficients fi and f2- When 
inserting the velocity potentials Eqs. (16a) and (17) into 
Eq. (14) for F'^^'^, we obtain a sum of terms each pro- 
portional to the square of 0^ = (j)^^ + (f)^^. This results in 
three types of contributions, (i) squares of (p^^ containing 
no information about the scattering and therefore yield- 
ing zero, (ii) squares of (p^^ proportional to the square of 
the particle volume and therefore negligible compared 
to {Hi) the mixed products (pij^(t>sc proportional to parti- 
cle volume a^, which are the dominant contributions to 
prad Keeping only these mixed terms representing in- 
terference between the incoming and the scattered wave, 
and using the index notation (including summation of re- 
peated indices), the ith component of Eq. (14) becomes 



^rad 



da Hi 



f (PinPsc)-Po(«^) 



I 

Jn 



dr d. 



(18a) 



f (Pi„Psc)-Po«"0 



I / in sc \ I / sc in \ 



(18b) 



=-/d.H 

Jn Po 



Po 



(Pin^iPsc) + (PsAPin) 

/ in o sc \ I / sc o in \ 



(18c) 



dr- - (Pindt^) - (PsAvT) 



(18d) 



= - /^dr l^{vrd,p,,)+Po{vrdjvf)^ (18e) 



(18f) 



Here, showing details not explained in Ref. 3, we have 

used = Cq Pi in Eq. (18a), Gauss's theorem in 
Eq. (18b), exchange of indices d^Vk = d^dj^cj) = df.d^4> = 
d^Vi to cancel terms in Eq. (18c), introduction of time 
derivatives by the continuity equation d^p^ = —PodjVij 
and the Navier-Stokes equation p^d^vi^i = —d^Pi = 
—CQdiPi in Eq. (18d), vanishing of time- averages of total 
time derivatives ((?((pwj)) = or (pQ^Wj) = — (f^Q^p) for 
cancelation and rearrangement in Eq. (ISe), and finally 
reintroduction of the vector potential (j)^^ in Eq. (18e). 

The d'Alembert wave operator — {l/c^)df acting 
on appears in the integrand of Eq. (18f), and since 



4>g^ is a sum of simple monopole and dipole terms, sig- 
nificant simplifications are possible. Just as the Laplace 
operator acting on the monopole potential (p = q/ (iircQr) 
yields the point-charge distribution, d^4> = ~{l/^o)^(''')^ 
in the static case, the d'Alembert operator acting on the 
retarded-time monopole and dipole expressions (17) also 
yields delta function distributions, 



1 



(19) 



= fi^d,p,^6{r) + h27ra'V- 



6{r) 



r » A. 



Now we see the great advantage of working in the far- 
field limit. The first term is easily integrated, when 
appearing in Eq. (18f), but for the second term we 
need to get rid of the divergence operator acting on 
the delta function before we can evaluate the integral. 
This we manage by Gauss's theorem. First we note that 

V • [?7(r)M(r)] = v'V ■ u + u ■ Vw for any scalar function 

V and vector function u. Therefore, Jg^ da n ■ (vu) = 
/j^ dr V • (vu) = dr {v'V ■ u + u ■ Vw), and we have 
derived the expression dr v'V ■ u — — J^dr u ■ Vv + 
Jg^da n ■ {vu). Now, since u oc v5{r) we obtain in 
Eq. (ISf) a volume integral encompassing the delta func- 
tion, thus yielding a non-zero contribution, and a surface 
integral avoiding the delta function, thus yielding zero. 
Consequently, the resulting expression for F'^'^ becomes 

A-TT 

(20a) 



(20b) 



47r 



3poCo 



2 (/iPinVPin) + 27ra3 Po{f2Vin-Vv,^) 



(20c) 



= —na 



2k, 



— ^ Re [/iP*,^Pin] - Po R-e [/2 < • VWi, 



with Pjj^ and v^^ evaluated at r = 0, (20d) 

where we have integrated over the delta function in 

Eq. (20a), applied the previously used rule (PiAt'^in) = 
-(^'in^tPin) in Eq. (20b), inserted p-^ = p-Jc^ and 
dfV^^ = — V^jj^/pq in Eq. (20c), and finally taken the 
time average using Eq. (11) in Eq. (20d). 

For standing waves (p^^, the spatial part /(r) of the 
incoming fields /(r)e'"* is real, so the nabla operator in 
Eq. (20d) does not lead to any phase changes. Conse- 
quently, the radiation force acting on a small particle 
(a ^ A) placed in a standing wave is a gradient force of 
the following form. 



^rad 



y^rad^ for a standing wave 



4n 



(21a) 
• (21b) 
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The radiation potential [/''^'^ is proportional to the vol- 
ume of the particle, and it contains a positive contribu- 
tion from the acoustic pressure fluctuations and a nega- 
tive contribution originating from the Bernoulli effect of 
the acoustic flow speed squared. 

For traveling waves, the spatial part /(r) of the incom- 
ing fields /(r)e"^* contains a phase changing factor, e.g. 
the plane-wave factor e''" '' or the spherical-wave factor 
e''^'", which changes the overall structure of the result- 
ing radiation force. Assuming for simplicity an incoming 
plane wave with wavenumber k parallel to v^^, we have 
Vpjjj = ifePin and = ikv, and Eq. (20d) leads to a 
resulting radiation force of the form 



^rad 



k, (22) 



for a purely traveling wave 
Note that this is not a gradient force. 



its monopole and dipole term becomes 
</.,,(r, 6) = </.^p(r) + ^^^{r, 6), a + h5<r « A, (24a) 

<t>^^{r) = -h-^dtp.J-, (24b) 



3/3o 



cos^ 



(24c) 



Finally, outermost in the far-field region r ^ A, the fluid 
is inviscid and compressible with a time-retarded (j)^^. 

In first-order scattering theory the monopole and 
dipole parts of the problem do not mix: is the co- 
efficient in the monopole scattering potential 4>^^ from a 
stationary sphere in the incoming density wave p^^, while 
/2 is the coefficient in the dipole scattering potential (f)^^ 
from an incompressible sphere moving with velocity 
in the incoming velocity wave v-^^. 



A. The monopole scattering coefficient /i 



IV. THE SCATTERING COEFFICIENTS 

The scattering coefficients and of Eq. (17), are 
found by matching the pressure Pi and velocity of 
the fluid with the boundary conditions at the particle 
moving with the instantaneous velocity v^. In the fol- 
lowing we use a spherical coordinate system with unit 
vectors (e^,eg,e^) located at the instantaneous center 
of the particle. Due to the azimuthal symmetry of the 
problem, all fields depend only on r and 6, and the veloc- 
ities has no azimuthal component, v = v^e^ + Vgeg. The 
polar axis points along the instantaneous direction of 
the incoming velocity v-^^^, such that v-^^^ = v-^^^e^. By the 
azimuthal symmetry of the problem, the particle must 
also move in that direction, v,^ = v„e^, 



smU v-e 



(23a) 
(23b) 



As sketched in Fig. 1(b), the response of the fluid is 
different in three regions of space [25]. Just outside the 
sphere, in the so-called acoustic boundary layer given by 
a < r < a -|- 5(5, viscosity is important due to the in- 
creased shear gradients in the velocity fields, as discussed 
after Eq. (7). Moreover, the fluid appears incompressible, 
since the time it takes an acoustic wave to propagate 
across the boundary layer around the particle is much 
less than the oscillation period, (a -|- 5i5)/co <C l/uj or 
a + 5S ^ A. The first-order pressure and velocity fields in 
the viscous and incompressible acoustic boundary layer 
are denoted p^^ and v^^^, respectively. In the next region, 
the so-called near-field region with a + 55 < r <^ X, the 
fluid is inviscid and compressible, but (j)^^ depends on the 
instantaneous argument t and not the time-retarded ar- 
gument t — r/cQ. The scattering potential Eq. (17) with 



The presence of a stationary, compressible particle 
causes a mass rate dtm of fluid to be ejected, that would 
otherwise have entered the particle volume. To flrst or- 
der, the ejection is determined by the mass flux pg'^sc 
carried by the scattered wave through a surface d^l en- 
compassing the particle in the near-field region. For a 
spherical surface with surface vector n = e^, we obtain 



[ in 

JdQ "5 



(25) 



The factor 1/3 was introduced in Eq. (17) to make the 
particle volume Vp = (47r/3)a'^ appear here. The rate of 
ejected mass can also be written in terms of the rate of 
change of the incoming density po + Pi„ multiplied by the 
particle volume Vp as dtm = dt [(po + Pin)^] • Expressing 
this through the compressibility k = ~{l/V){dV/dp) of 
the fiuid, Kq, and the particle, ftp, we obtain. 



Vp dcp^^ 



(26) 



Now, fi is obtained by equating Eqs. (25) and (26), 



(27) 



This result is identical to that of Gorkov [3] . We note that 
f-^ is real- valued and depends only on the compressibility 
ratio k between the particle and the fluid; the viscosity 
of the fluid does not influence the compressibility and 
mass ejection. For identical compressibilities, k = 1, the 
monopole scattering vanishes, /i(l) = 0. 



B. The dipole scattering coefficient /2 

As /2 is related to the translational motion of the 
particle, it depends on the viscosity of the fluid, and 
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we must therefore explicitly calculate the first-order ve- 
locity v^^(r,9) in the viscous, acoustic boundary layer 
a < r < a -f- 55. This velocity must match the dipole 
part Vjj^ -I- V(/)^p of the fluid velocity in the near-field re- 
gion r ^ A, see Eqs. (23a) and (24c). Because of the 
separation of length scales the matching can be made at 
a distance r « r* fulfilling a + 5(5 <C r* ^ A, and this 
so-called asymptotic matching [25, 31] is written 



2-'2 ^3 



Vi„cos6», (28a) 
v,^{-sm0). (28b) 



At the surface of the sphere, r = a, the no-slip boundary 
condition requires v^^ to equal Vp of Eq. (23b) , 



^ab(a,fi') = Vpcose, 
^ab(a,^') = Wp(-sin( 



(29a) 
(29b) 



The velocity of the sphere is given by Newton's second 
law with dfVp = — iwUp and with the viscous stress from 



the fluid acting on the surface of the sphere, 



- IgTT a PpLO Vp = 



dan - (T. 



ab 



(30) 



= 2na^ J d(c 



cos 0) ( - + Kr) COS e - sin 9 



The viscous stress components are a^J? = 2r]d^vf^ and 
< = ^[(l/r)a,«f + drvf - {l/r)vf] . 

The determination of the pressure p^^^ and velocity 
field is eased by the incompressibility of the fluid 
in the acoustic boundary layer, V • v^y^ = 0. Tak- 
ing the divergence of the first-order Navier Stokes equa- 
tion (6b), leads to the Laplace equation of the pressure, 
V^Pjjj, = 0. Since we are seeking the dipole solution, and 
since p^y^ must match asymptotically with the dipole part 
ipow((/>in -|- 0jp) of Eqs. (16c) and (24), we can immedi- 
ately write down the pressure in the boundary layer. 



la3 

2 r2 ^2 



V- cos f 



Pe.hi"'^^) =^P0^a\l + -f2 



1 



V;„ COS 6. 



(31a) 
(31b) 



For the velocity field, incompressibility combined with 

the azimuthal symmetry of the problem, implies that v^^ 
can be written in terms of a stream function ^'(r, 9) as 



r sm y r 



(32a) 
(32b) 



Taking the rotation of the first-order Navier-Stokes equa- 
tion (6b) and using Eq. (32a) leads to an equation for ^, 



(V^ + g2) [^^,^r, 9) e J = 0, with q : 



(33) 



The solution to this biharmonic vector equation can be 
found as the sum 'i/ = '^^ + 'i>2, where 

V2(*ie^) = V^*! - = 0, (34a) 

^ sin 9 

*2 2t 
2 2n = -«^2- 

sm 9 

(34b) 



The dipole part of the solution to Eq. (34a) is the Leg- 

endrc form '^i{r,9) = A^r cos + A2 cos 9 / r"^ . Like- 
wise, the dipole part of the solution to Eq. (34b) is 
the Hankel form ^'2(7',^) = Bh\{qr) auj^^sin^, where 
^li^) ~ — e'*(s + i)/s^ is the spherical Hankel function 
of the first kind (outgoing wave) of order 1. The con- 
stants A-y and A2 are given by the asymptotic matching 
conditions (28), so ^ and v^^ = V x (^'e ,) become 



lr-^'^ + ahl{qr)B 
l-/24 + 2,aB^^^(-^) 



V- sm 



2 r3 



1 + ^— + qaB 



(35a) 
Win cos 61, (35b) 

?;;„(- sin 61). 
(35c) 



Note that all information about the viscous boundary 
layer is given through the constant q = {1 + i)/6 of 
Eq. (33). Since the spherical Hankel function decays ex- 
ponentially on the length scale S, hl{qr) oc e'*'' oc e~''^, 
the viscous, acoustic boundary layer has a width of ~ 56. 
An important, dimensionless parameter is the ratio S of 
the viscous penetration length S and the particle radius a, 



-6='-. 

a 



(36) 



Inserting v^^^ of Eqs. (35b) and (35c) into the no-slip 
condition at the surface of the sphere, Eq. (29), and the 
pressure ^Jg^jj of Eq. (31b) together with v^y^ into Newton's 
second law for the sphere, Eq. (30), we arrive at three 
equations for the three unknowns /2, Vp, and B, 



l-f2 + 2Bhliqa) 



l + -f^ + Bd,{shl{s)) 
l + lf2 + iS^Giqa)B 



qa 



(37a) 
(37b) 

(37c) 



where the auxiliary function G{s) is given by 



G{s) = 2s'd, 



s 



2h\ 



sd^ {sh\ )-2h\= is(i + s) hi (s). 



-dAsh\) 
(38) 



with hl{s) = -(i/s)e'''. When subtracting Eq. (37a) 
from Eq. (37b), d^{sh\{s)) - 2h\{s) = shl{s) - 3hl{s) 



appears. It is therefore useful to introduce the dimen- 
sionless, i5-dependent variable 7 given by 



7(^) 



3h\ (qa) 



l + i(l + ^) 



d. 



qahliqa) 

By straightforward algebra, /2 is found from Eq. (37), 
2[l-7(^)](p-l 



(39) 



2p+l-37((5) 



(40) 



The viscosity-dependent dipole scattering coefficient 
is in general a complex- valued number, and its real and 
imaginary values are abbreviated as 



f^{p,S)^lm [f,{p,~S)]. 



(41a) 
(41b) 



In the absence of viscosity, S 
valued result by Gorkov [3], 



0, we recover the real- 



/2(P,0) 



2(P-1) 
2p+l ' 



(42) 



Physically, this result for an inviscid fluid can also be de- 
rived directly from Eq. (37); The acoustic boundary layer 
vanishes, B = 0, and the condition (37b) on the tangen- 
tial velocity component is dropped, so we are left with 
the normal-component condition (37a), = {1 — f2)''^in^ 
and Newton's second law (37c), pv^ = (l + ^./2)^in' from 
which Eq. (42) follows. 

The non-zero imaginary part of for 6 > implies 
that for a traveling wave (f)^j^ the product terms (/'i„0sc of 
order remain finite. This is in contrast to the inviscid 
case, where these terms vanish and the radiation force 
is reduced by a factor of (ka)^ — (27ra/A)^ since only 
the quadratic terms 0^^, proportional to remain finite 
[2, 3]. In agreement with Doinikov [26], our analysis thus 
predicts the possibility of realizing a radiation force for 
traveling waves, which is a factor of (fca)"'^ — (A/27ra)^ 
stronger than expected from the standard inviscid theory. 



C. Properties of 

As the dipole scattering coefficient /j, in contrast to 
the monopole coefficient fi, depends on viscosity, we 
study some of its properties in more detail below, see 
Fig. 2. Importantly, /2 is zero for neutral-buoyancy par- 
ticles {p = 1) irrespective of the viscosity. 



/2(1,^)=0, 



(43) 



and generally small for near-neutral-buoyancy particles. 

For a large particle in a low- viscosity fluid S 1, the 
correction to the Gorkov expression for is found by 
Taylor-expanding Eq. (40) to flrst order in 5, 



/2(P,<5«1) 



2(P-1) 
2p+l 



, (44) 



[a) 


1 K 
1. /J 




1.50 - 


0" 


1.25 ■ 








1.00 1 




0.75 - 








0.50 ■ 




0.25 ■ 



0.00 





p = 2.00 ■ 




p = 1.75 


p = 1.50 




p = 1.25 ■ 


j^-""' p = 1.00 


^V~^ 5 - 0.75 




p = 0.50 - 


— - 0.25 




3 = 0.00 






FIG. 2. (a) The real part of the viscosity-dependent 
dipole scattering coefficient relative to its inviscid counter- 
part f2 {p, S) / f2 (p,0) plotted versus the non-dimensionalized 
thickness S of the viscous, acoustic boundary layer, (b) The 
imaginary part f2ip,S) versus 5 for p — 1.25, 1.50, 1.75, and 
2.00. For the same parameters 0.17 < f^ip, 00) < 0.67. 



in agreement with Doinikov [26]. In earlier work by 
Weiser and Apfel [32] (building on Urick [33]) the nu- 
merator of the viscous correction was found to be (9/2)5 
instead of 3{p — 1)6. However, this discrepancy is due to 
an imprecise treatment of the viscous boundary layer in 
the earlier work. 

For a small particle in a high- viscosity fluid S ^ 1, the 
viscosity dependence saturates, as the particle essentially 
becomes a point singularity at the center of the acoustic 
boundary, and becomes 



/2(/5,(5> 1) « - (/5- 



!)■ 



(45) 



Doinikov presented the same expression except for the 
overall sign. This may be a misprint in his paper, as the 
consequence would otherwise be an un-physical reversal 
of the sign of the force as S is increased from zero to 
infinity. 

In Fig. 2(a) is shown plots of the real part /2 (p, S) of 
the viscosity-dependent dipole scattering coefficient rela- 
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tivc to the inviscid coefficient f2 {p, 0) as a function of i5 
for different values of p between and 2. For these values 
of p the values of f2{p, <5)//2^(Pi 0) fall in the range from 
0.3 to 1.7, and the saturation of f2 sets in for moderate 
values of 6 between 1 and 2. 

In Fig. 2(b) is shown the imaginary part f2{p,^), 

f ^^ ^ 6{i-pni+s)~s 

' (l + 2p)2+9(l + 2p)5+f (^2 + ^3+ 1^4)' 

(46) 

with 

6(1 - p)2 - 
/.■(P,^«l)^ (;^^g, 6, (47a) 

_ 94 _ 

f^{p,6^1)^-{l-pr6-\ (47b) 

It exhibits a marked maximum for S w 0.5 with an ampli- 
tude roughly one order of magnitude smaller than the sat- 
uration value of f2{p, oo), Eq. (45), for the corresponding 
densities p. 



TABLE I. List of values of material parameters at 20 °C for 
typical liquids [water (wa), NaCl solution (scs), percoU (pc), 
glycerol (gl)] and solids [pyrcx (PY), polystyrene (PS), poly- 
methaerylatc (PM), melaminc resin (MR), a representative 
biological cell (Cell)] used in microchannel acoustophoresis. 



Material 


Density 


Compress. 


Longitud. speed 


Viscosity 




P [kg/m^] 


K [1/TPa] 


of sound c [m/s] 


r] [mPa-s] 


wa'' 


998.2 


456 


1482 


1.002 


scs'' 


1071 


365 


1599 


1.170 


pc^ 


11.30 


390 


1507 


100 




i2()i 


219 


1901 


1112 




2230 


27.8 


5674 




PS'= 


1050 


172 


2350 




PM'^ 


1190 


148 


2380 




MR'^ 


1510 


67.5 


3132 




Cell"* 


1100 


400 


1500 





^ Prom Ref. 35 

Sodium chloride solution of salinity S = 0.1, see Appendix A 
" Prom Sigma-Aldrich Production GmbH and Fluka data sheets 
<* Prom Refs. 14 and 36 



D. Resulting expressions for the radiation force 

In summary, the main result of the paper are the fol- 
lowing analytical expressions for the acoustic radiation 
force F'^'^ on a spherical particle of radius a, density Pp, 
and compressibility Kp suspended in a fluid of density 
Pq, compressibility Kq, and viscosity rj, and exposed to a 
first-order standing and traveling acoustic wave p^^ and 
Ujj^ in the long wavelength limit A a. 

For a standing acoustic wave we have obtained 



jjrad 



rad 



47r 



1 3 



7(5) 



1 — K, with K = — , 



(48a) 
(48b) 

(48c) 



= Re 



2[l-7(^)](/5-l) 



, with p = 



Po' 



6, with 6 = —, 
a 



(48d) 

(48e) 



2p + 1 - 37(5) 
1 + i(l + 6) 

and for a traveling planar wave with wavevector fe, 

F^^ = f\{p,l)T,c?p^{v1^)k, (49) 
where /2 (p, 5) is given by Eq. (46). 

V. EXPERIMENTAL IMPLICATIONS 

A. Typical materials 

In practical applications, especially involving biologi- 
cal samples, the solvents are often aqueous salt solutions. 



Among these, sodium chloride (NaCl) solutions are ar- 
guably the ones best characterized acoustically, so we use 
this solvent as one of our model liquids in the following 
analysis. In Appendix A the current best values of the 
speed of sound Cg^.^, the density p^^g, and the viscosity t]^^^ 
of sodium chloride solutions (scs) are given as a function 
of temperature T (in °C) and mass fraction S of NaCl in 
the solution (salinity). To study the effects of change in 
viscosity we also include calculations with glycerol (gl) 
and with pcrcoU (pc), a solution of polyvinylpyrrolidone- 
coated silica nanoparticles for which the speed of sound 
is nearly the same as for pure water [34], see Table I. 

As typical materials for the particles we have chosen to 
analyze the polymers polystyrene (PS), polymethacrylate 
(PM), and melamine resin (MR), as well as pyrex (PY) 
and a typical biological cell (Cell), see Table I. 



B. Traveling acoustic waves 

As mentioned above, our theory including viscosity 

predicts a strong enhancement of the acoustic radiation 
force by a factor of (ka)^^ compared to the standard 
inviscid theory for the case of purely traveling waves. 
Here, we analyze the case of a planar traveling wave with 
fe = fcCj, and = Pjje'('^^~"*\ In this case the acoustic 
energy density E^^ 
radiation force (22) becomes 



of the wave is E^^ = ^KqpI, and the 



nrad 



= -KO?kf^{p, S)E^e^ 



(50) 



For a 5-]im-diameter pyrex sphere in water at 1.5 MHz 

wc have a = 2.5 pm, S = 0.23, k = 4.24 x 10^ m-\ 
and p = 2.23. From Eq. (46) we find /J (2.23, 0.23) = 
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0.058, and taking the typical acoustic energy density 
E^^ = 100 J/m^ [22], we arrive at F""^ « 1.2 pN. Under 
the influence of this force, the pyrex sphere would reach 
the terminal translational velocity = F'^^'^ / (Qnrja) m 
25 jim/s. This is a significant velocity for microchannel 
acoustophoresis, where typical velocities lie in the range 
from 5 to 500 pm/s [23]. In the standard inviscid theory 
[2, 3], the estimate for the radiation force is a factor of 
{kaY ~ 10~^ lower, corresponding to an acoustophoretic 
velocity smaller than 0.1 nm/s. 



C. Standing acoustic waves 

In many experiments on rectangular microfluidic chan- 
nels with coplanar walls at z = and z = h, the in- 
coming wave have approximately been a resonant, stand- 
ing ID pressure wave of the form pi = p^cos{kz), with 
wavenumber k — mr/h, where n is the number of half 
wavelengths, and with the acoustic energy density E^^ = 
■jKqP^. The expression for the radiation force then sim- 
plifies to the classic result by Yosioka and Kawasima [2] , 



F: 



rad 
ID 



47r $(«;, p, (5) a^fcE'a^j. sin(2fcz), (51a) 



(51b) 



where the acoustophoretic contrast factor $(k, p, (5) now 
depends on viscosity. 

Experiments on suspended biological cells involve near- 
ncutral-buoyant particles, \p— 1| ^ 1, implying that the 
monopole coefficient \fi\ is typically much larger than 
the dipole coefficient 1/2"^!. Because the acoustic contrast 
factor $ defined in Eq. (51b) is a linear combination of 
the two scattering coefficients, a good quantitative mea- 
sure of the ability to detect the effect of viscosity on the 
acoustic radiation force is therefore the relative change in 
$ with and without viscosity. We therefore find it help- 
ful to introduce the detectability measure T) of viscous 
effects as 



~ ?x ^)-$(k,p,0) 



$(k,p,0) 



(52a) 
(52b) 



Examples of 2? for in NaCl solutions are shown in 
Fig. 3 going from nearly undetectable sub-l-%-levels for 
polystrene spheres to above 10-%-levels for pyrex spheres. 

The effect of including the viscosity in the expression 
for the acoustic radiation force can also be illustrated by 
contour plots of the contrast factor $ in the (k, p)-plane 
for fixed values of i5 as shown in Fig. 4. The change in the 
contrast factor is clearly seen by the changing contour 
lines. While $ is independent of S along the neutral- 
buoyancy line p = 1, its value is increased when going 
from the inviscid case 5 = in panel (a) to the viscous 
case (5 = 1 in panel (b) . The change of the contour line 




5 10 15 20 25 

NaCl concentration (mass fraction S in %) 

FIG. 3. The detectability V for solid microspheres of radius 
a ranging from 0.1 11m to 10 ym in a sodium chloride solution 
as a function of NaCl concentration (mass fraction S), see 
Appendix A. (a) Polystyrene spheres, for which 2? is typically 
a few percent or lower, (b) Pyrex spheres, for which D easily 
can be larger than 10%. 



$ ~ 0.0 is particularly interesting as particles on op- 
posite side of this line move in opposite directions, and 
the plot of $ in the (k, /5)-plane is therefore also useful 
when attempting to tune the solvent to obtain binary 
separation of particles. In Fig. 4, a number of specific 
examples of materials are marked by crosses. As particle 
material are chosen polystyrene (PS), polymethacrylate 
(PM), melamine resin (MR), and a typical biological cell 
(Cell), while the liquids are water (wa), glycerol (gl), and 
percoU (pc). Note that the Cell/gl and Cell/wa points lie 
on opposite sides of the zero contour. A curve connecting 
these two points would represent the acoustophoretic re- 
sponse of cells in various mixtures of glycerol and water. 
From a purely physical point of view, this system may 
therefore form an excellent tunable solvent with respect 
to obtaining binary separation of cells. 

The effect of viscosity can also be studied through 
vertical acoustic trapping [37, 38], where the buoy- 
ancy force F'^"°y = (47r/3)(p - l)pQa^g is balanced by 
a vertically-oriented standing plane-wave acoustic field 



10 




0.5 



0.0 



2.0 



1.5 



MR/wa 
PS/wa \ PY/gl 



PM/wa it 



PY/pc 



0.5 



. PY/wa 

X IJS 





MR/wa 

PS/wa \ PY/gl 

" \ \ PY/pc 
PM/wa ^ \ 



PY/wa 0.7 



1.0 



1.5 



2.0 



2.5 



FIG. 4. Contour plots of the acoustic contrast factor ^{k, p, 5) 
as function of k and p for fixed values of 5, from <1> < —0.5 
(black) to $ > 0.5 (white) in steps of 0.1. The position in 
the (k, p)-plane for various material parameters are marked 
by crosses and the following labels: polystyrene (PS), poly- 
methacrylate (PM), melamine resin (MR), atypical biological 
cell (Cell), water (wa), glycerol (gl), and percoU (pc). (a) The 
inviscid case (5 = 0. (b) The viscous case with 5 = 1. 



prad _ ^radg^^ p^j. ^ given acoustic energy density E^^, 
the maximal acoustic radiation force is given by the am- 
plitude A7r^{K, p, 6)a^kE^^ in Eq. (51a). Based on this, 
the critical trapping force is defined as the threshold for 
obtaining vertical acoustic trapping using the smallest 
possible acoustic energy density i?™'". 



II 



Po9 



3$(k,p,(5) k 
The effect of viscosity can therefore be measured as 

1 



(53) 



l + V 



(54) 



A quantity more readily accessible experimentally may 
be the voltage 11^^ used to drive the piezo transducer 
generating the ultrasound wave in a typical experiment. 
As E^^ scales with the square of U^^ [22], we have 



t/pz(0). 



VTTv 

Here, the detectability V of Fig. 3 appears directly. 



(55) 



VI. CONCLUSION 

We have derived an analytical expression for the acous- 
tic radiation force in the long- wavelength limit (5, a <C A 
on a compressible, spherical particle of radius a sus- 
pended in a liquid with viscous penetration depth S. 

We have analyzed the experimental predictions pro- 
vided by our expression for traveling waves and for stand- 
ing waves. In the case of the former we find a strong en- 
hancement proportional to {ka)~^ ~ 10^ relative to the 
inviscid case due to non- vanishing interference between 
the incoming wave and the scattered wave. For standing 
waves, we have found a negligible sub-1% deviation from 
the inviscid result for large (micrometer-sized), nearly- 
neutral-buoyancy particles, such as biological cells in wa- 
ter. However, significant deviations above 10% from the 
inviscid result were found for buoyant particles, such as 
pyrex in water. The smaller the particle radius a is rela- 
tive to the viscous boundary-layer thickness S, the larger 
the effect of viscosity. 

It should be possible using state-of-the-art instrumen- 
tation for microchannel acoustophoresis, such as the au- 
tomated micro-PIV setup recently published by Augusts- 
son et al. [23] , to experimentally test the predictions pre- 
sented in this paper. 
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Appendix A: Material parameters of NaCl solutions 

In the following we summarize parameter values for the 
speed of sound c^^^ , the density p^^^ , and the viscosity rj^^^ 
of sodium chloride solutions (NaCl in water) as a function 
of temperature T (in °C) and mass fraction S of NaCl in 
the solution (salinity). The salinity ranges from zero in 
pure water to a maximum of 0.26 in a saturated solution. 

The speed of sound c^csi^^ ^) NaCl/water solutions 
is given by Kleis and Sanchez [39] as 



4 

E 

J=0 



(a - +b,S) 



T 



1 °C 



(Al) 



11 



where the coefBcients and bj are 





1.40309 X 10'\ 


6o = 


1.40190 X 10^ (A2) 


ai = 


4.68391 X 10°, 


bi = 


-1.14996 X 10\ 


^2 = 


-4.05388 X 10-2, 


b,= 


2.23748 X 10"^, 


Og = 


1.29550 X 10-^ 


63 = 


1.48238 X 10"^, 


a4 = 


6.91485 X 10"^ 


&4 = 


-9.46165 X 10"^. 



The density p^^JSjT) of sodium chloride solutions is 
given by Laliberte and Cooper [40] as 



l-5 + 5Kpp(5,TKa(T)' 
where the density p^a^(T) of water is given by 



(A3) 



Pwa(r) = 



5 



1 + d-^ ^ 



1 



kg m 



(A4) 



with the coefficients d and dj being 



d = 1.6879850 x 10" 



dn 



9.9983952 x 10^ 



rfi = 1.6945176 X 10\ d^ = 
d2 = -7.9870401 x 10-^ 4 = 



(A5) 

-4.6170461 X 10"'\ 
1.0556302 X 10"^ 
-2.8054253 x 10"^°, 



and where the apparent specific volume V {S,T) of 



NaCl is given by 



I? + 62 + 63 



5(1^ + ^4) 



(A6) 



{cqS + e^) exp 

with the coefficients e^- being 

eo = -4.330 X 10-^ 63 = 1.4624 x 10-^ 
ei = 6.471 X 10-^ 64 = 3.3156 x 10^ 
62 = 1.01660 X 10°, eg = 1.0000 x 10"^ 

The viscosity ri^^^{S,T) of NaCl/ water solutions is 
given by Laliberte [41] as 

r?,,,(5, T) = [n^,{T)] ^'-'^ [,7NaCi(5, T)]', (A7) 

where the viscosities rj^^{T) and %aCi('^' water and 
liquid NaCl, respectively, are 

(J.. [t^ + ^4(>; mPa s 

'"''^ [0.05594 + 5.2842] + 137.37 



'?NaCi('S'.^) =exp 
with the coefficients being 



(A8) 

mPas, 
(A9) 



= 1.6220 X 10\ /i4 = 7.4691 x 10^^ 
/12 = 1.3229 X 10°, /i5 = 3.0780 x 10\ 
/13 = 1.4849 X 10°, Hq = 2.0583 x 10°. 
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